from deslab import *
from random import sample 
import time
from subprocess import check_call
import string
from exp_diagnosis import is_diagnosable_verif as test_verifier, nonobservable_scc
from exp_diagnosis import is_diagnosable_diag as test_diagnoser,diagnoser
from collections import Counter
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from scipy import stats
from matplotlib import rc
import math
from numpy import linalg
import numpy as np
import networkx as nx
import shelve

#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)


INSTANCES = [(3,2),(3,3),(3,4),(3,5),(4,2),(4,3),(5,2)]
SIGMA = dict(zip(range(0,27), string.ascii_lowercase,))
syms('a b c d e f g')
SIGOBS = {2: set([a]), 
          3: set([a,b]),
          4: set([a,b,c]),
          5: set([a,b,c,d]),
          6: set([a,b,c,d,e]),
          7: set([a,b,c,d,e,f])
          }

COLOR={2:(0,0,0.75),3:(0,0.5,0),4:(0.5,0,0),5:(0,0,0.75)}

RESULTFILE = '../instances/resultexh_c'

"""DESCRIPTION

In this experiment we exhaustively generate all the transition structures with n states
and k events, the failure event is the last event in the alphabet being considered"""



def generate_instances(n,k):
    c_i = 0
    S = exh_incomplete_fsa(n,k)
    tt=time.time()
    G_S = []
    while S: 
        #G_S.append(S[:]) 
        S = exh_incomplete_fsa(n,k)  
        if c_i % 500000 ==0: 
            print c_i,  time.time()-tt   
        c_i +=1
    print "writing file to hard disk"
    save(G_S,'../instances/Gexhaustive-'+str(n)+'n-'+str(k)+'k.exp')
    save(c_i,'../instances/sizeofGexhaustive-'+str(n)+'n-'+str(k)+'k.exp')
    print "file transfered to disk "
    print c_i 
   
   
#def exaustivegen(len_X,len_S,num_unobs):
#    tt=time.time()
#    c_i = 0
#    len_V = []
#    len_D = []
#    test = []
#    S = [SIGMA[e] for e in range(len_S)]
#    X = range(len_X)
#    So = [SIGMA[e] for e in range(len_S-num_unobs-1)]
#    X0 = [0]
#    Xm = [] 
#    sf = SIGMA[len_S-1]
#
#
#    for T in Gv['transitions']:
#        G = fsa(X, S, T, X0, Xm, Sigobs = So) 
#        #print "calculating diagnoser"
#        Gd = diagnoser(G,sf)      
#        if Gd.name:            
#            len_D.append(len(Gd))
#        t, V = test_verifier(G,sf)       
#        len_V.append(len(V))
#        #print len(Gd)
#        test.append("W")
#        if c_i%10 ==0:
#            print c_i, time.time()-tt
#            print "mean: ",np.mean(len_D),"curr: ",len(Gd)
#        c_i +=1           
#    result_exp = {'n': n,
#                  'k': k,
#                  'instances':i,
#                  'valid': len(len_V),
#                  'len_D': len_D,
#                  'len_V': len_V,  
#                  'ind_G': ind_G,  
#                  'test' : test,
#                  'max_D': max(len_D),
#                  'max_V': max(len_V)                 
#                  }
#
#    res_db=shelve.open(INSTANCEFILE+"_1u.db")
#    print "loading results of experiment %sX-%sS-%sO"%(len_X,len_S,num_unobs)
#    res_db['%sX-%sS-%sO'%(len_X,len_S,num_unobs)] = result_exp
#    res_db.close()
#    print " results loaded" 
    
def entropy(G):
    """This function calculates the entropy of automaton G"""
    # the graph $$ is a living object
    adj_matrix = nx.to_numpy_matrix(G.Graph)
    eigval,eigvec = linalg.eig(adj_matrix)
    abs_eigval = abs(eigval)
    maximum = abs_eigval.max()
    entropy = maximum#math.log(maximum,2)    
    return entropy


def calculate_diagver(n, k,num_unobs):
    i = 0
    splitnum = 8000000
    len_V = []
    len_D = []
    ind_G = []
    sf = SIGMA[k-1]
    test=[]
    tt=time.time()
    c_s = 1
    S = exh_incomplete_fsa(n,k)
    So = [SIGMA[e] for e in range(k-num_unobs-1)]    
    
    while S:
        G = istr2fsa(S,n,k) 
        G = G.setpar(Sigobs=So)       
        nonlive = [x for x in G if G.Gamma(x) == set()]
        if i % 1000 == 0: print i, time.time()-tt  
        if not nonlive:
            if not nonobservable_scc(G):
                if any(sf in G.Gamma(x) for x in G):                
                    Gd = diagnoser(G,sf)                    
                    t, V = test_verifier(G,sf)
                    len_V.append(len(V))
                    len_D.append(len(Gd))
                    ind_G.append(i)               
                    test.append(t)
        S = exh_incomplete_fsa(n,k)  
        i +=1
        if i % splitnum==0:
        
            result_exp = {'n': n,
                          'k': k,
                          'instances':i,
                          'valid': len(len_V),
                          'len_D': len_D,
                          'len_V': len_V,  
                          'ind_G': ind_G,  
                          'test' : test,             
                          }
            res_db=shelve.open(RESULTFILE)
            res_db['%sn-%sk-%su-%sp'%(n,k,num_unobs,c_s)]=result_exp
            print "loaded instances in resultexh[%sn-%sk-%su-%sp]"%(n,k,num_unobs,c_s)
            res_db.close()
            len_V = []
            len_D = []
            ind_G = []
            test = []
            c_s +=1
            
            # nada que actualiza esta tortas            
    result_exp = {'n': n,
                  'k': k,
                  'instances':i,
                  'valid': len(len_V),
                  'len_D': len_D,
                  'len_V': len_V,  
                  'ind_G': ind_G,  
                  'test' : test,         
                  }
    res_db=shelve.open('../instances/resultexh.db')
    res_db['%sn-%sk-%su-%sp'%(n,k,num_unobs,c_s)]=result_exp
    print "loaded instances in resultexh[%sn-%sk-%su-%sp]"%(n,k,num_unobs,c_s)
    res_db.close()

  
def generate_graphics(color, len_X, len_S, viz =False):    
    R = load('../instances/result-'+str(len_X)+'n-'+str(len_S)+'k.exp')
    print "result loaded"
    max_D = R['max_D']
    max_V = R['max_V']
    ind_max_D = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==max_D]
    ind_max_V = [R['ind_G'][j] for j, vm in enumerate(R['len_D']) if vm==max_V]
    num_max_D = len(ind_max_D)
    num_max_V = len(ind_max_V)
    num_diag = len([i for i in R['test'] if i==True])        
    
    print 'valid', len(R['ind_G'])
    print 'diagnosable', num_diag
    print 'max D', max_D
    print 'max V', max_V
    print 'mean D',  np.mean(R['len_D'])
    print 'mean V',  np.mean(R['len_V'])
    print 'maximum worst case D', len(ind_max_D)
    print 'maximum worst case V', len(ind_max_V)
    print "$(%s,%s)$ & %s & %s & %s & %s & %s & %.2f  & %.2f \\\\ "%(len_X, len_S, R['instances'], len(R['ind_G']), num_diag, max_D, max_V, np.mean(R['len_D']), np.mean(R['len_V']))
    
    H = Counter(R['len_D'])
    x = [i for i,j in H.items()]
    y = [j for i,j in H.items()]
    
    #V = Counter(R['len_V'])
    #xv =  [i for i,j in V.items()]
    #yv =  [j for i,j in V.items()]
    # There are n+1 possible number of "successes": 0 to n.
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)
    #l2=plt.step(xv,yv,".",linewidth=1.25,color=(0.5,0.5,0.5,0.5))
    plt.plot(x,y,"bo", color=color)    
    l1=plt.vlines(x, [0], y, color=color,linewidth=1.5)
    # Format x-axis and y-axis.
    plt.axis([0, max(x)*1.05, 0, max(y)*1.05])
    plt.xticks(range(0,max_D+1,5))
    plt.title(r"Frequency of state complexity of $G_d$ for $||X||=%s$ and $||\Sigma||=%s$"%(len_X,len_S))
    plt.xlabel("Number of states of diagnoser\\\\ (worst theoretical case $2^{2\\times%s}=%s$)"%(len_X,2**(2*len_X)))
    plt.ylabel("Valid instances")
    #plt.legend([l1],["$G_d$"])
    f_path='../document/figures/Figure-freq-%sn-%sk.pdf'%(len_X, len_S)
    plt.savefig(f_path, dpi=None, format='pdf')
    print "loaded Figure-freq-%sn-%sk.pdf"%(len_X, len_S)
    if viz:
        command = 'evince'+' '+f_path+'&'
        retcode = check_call(command, shell=True, stderr=0)



def generate_graphics_diagnosable(color): 
    R = load('../instances/result-'+str(n)+'n-'+str(k)+'k.exp')
    print "result loaded"
    max_D = R['max_D']
    
    
    
    ind_max_D = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==max_D]
    #ind_max_V = [R['ind_G'][j] for j, vm in enumerate(R['len_D']) if vm==max_V]
    #num_max_D = len(ind_max_D)
    #num_max_D = len(ind_max_V)
    ind_diag = [i for i,j in enumerate(R['test']) if j==True]
    
    D_I = [R['len_D'][i] for i in ind_diag]
    
    #print 'total', len(S)
    #print 'valid', len(R['ind_G'])
    #print 'diagnosable', num_diag
    #print 'max D', max_D
    #print 'max V', max_V
    #print 'mean D',  np.mean(R['len_D'])
    #print 'mean D',  np.mean(R['len_V'])
    #print 'maximum worst case D', len(ind_max_D)
    #print 'maximum worst case V', len(ind_max_V)
    #print "$(%s,%s)$ & %s & %s & %s & %s & %s & %.2f  & %.2f \\\\ "%(n,k,len(S), len(R['ind_G']), num_diag, max_D, max_V, np.mean(R['len_D']), np.mean(R['len_V']))
    
    
    
    H = Counter(D_I)
    x = [i for i,j in H.items()]
    y = [float(j)/len(D_I) for i,j in H.items()]
    
    
    # There are n+1 possible number of "successes": 0 to n.
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)
    plt.plot(x,y,"bo", color="blue")
    plt.vlines(x, [0], y, color=(0,0,1), linewidth=1.5)
    # Format x-axis and y-axis.
    plt.axis([0, max(x)*1.05, 0, max(y)*1.05])
    plt.xticks(x)
    plt.title(r"Probability of state complexity of diagnoser for $n=%s$ and $k=%s$"%(n,k))
    plt.xlabel("Number of states of diagnoser\\\\ (worst theoretical case $2^{2\\times%s}=%s$)"%(n,2**(2*n)))
    plt.ylabel("Probability")
    plt.savefig('../document/figures/Figure-diagnosable-%sn-%sk.pdf'%(n,k), dpi=None, format='pdf')
    
    #plt.show()



def analize_diagnosability(color, len_X, len_S, viz =False):    
    R = load('../instances/result-'+str(len_X)+'n-'+str(len_S)+'k.exp')
    print "result loaded"

    max_D = R['max_D']
       
    H = Counter(R['len_D'])
    x = [i for i,j in H.items()]
    y = [j for i,j in H.items()]
    y_diag =[]
    x_diag=[]
    for ci in x:
        ind_si = [i for i, dm in enumerate(R['len_D']) if dm==ci]        
        si_diag = len([i for i in ind_si if R['test'][i]==True])
        if si_diag > 0:
            x_diag.append(ci)
            y_diag.append(si_diag)
                
    color2=(0.8, 0, 0.8)    
    #V = Counter(R['len_V'])
    #xv =  [i for i,j in V.items()]
    #yv =  [j for i,j in V.items()]
    # There are n+1 possible number of "successes": 0 to n.
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)
    #l2=plt.step(xv,yv,".",linewidth=1.25,color=(0.5,0.5,0.5,0.5))
    plt.plot(x,y,"bo", color=color)    
    l1 = plt.vlines(x, [0], y, color=color,linewidth=1.5)
    l2 = plt.plot(x_diag,y_diag,"_", color=color2,linewidth=1.5,markersize=10)   
    l3 = plt.vlines(x_diag, [0], y_diag, color=color2,linewidth=1.5)
    # Format x-axis and y-axis.
    plt.axis([0, max(x)*1.05, 0, max(y)*1.05])
    plt.xticks(range(0,max_D+1,5))
    plt.title(r"Frequency of state complexity of $G_d$ for $||X||=%s$ and $||\Sigma||=%s$"%(len_X,len_S))
    plt.xlabel("Number of states of diagnoser\\\\ (worst theoretical case $2^{2\\times%s}=%s$)"%(len_X,2**(2*len_X)))
    plt.ylabel("Valid instances")
    plt.legend([l3],["Diagnosable instances"])
    f_path='../document/figures/Fig-diagnosiscomp-%sn-%sk.pdf'%(len_X, len_S)
    plt.savefig(f_path, dpi=None, format='pdf')
    print "loaded Figure-diagnosiscomp-%sn-%sk.pdf"%(len_X, len_S)
    if viz:
        command = 'evince'+' '+f_path+'&'
        retcode = check_call(command, shell=True, stderr=0)
        




def analize_entropy(color, len_X, len_S, viz =False):    
    R = load('../instances/result-'+str(len_X)+'n-'+str(len_S)+'k.exp')
    print "result loaded"
    S = load('../instances/Gexhaustive-'+str(len_X)+'n-'+str(len_S)+'k.exp')
    print "instances result loaded"
    max_D = R['max_D']         
    C = Counter(R['len_D'])
    x = [i for i,j in C.items()]
    ymax = max([j for i,j in C.items()])
    y = [float(j)/ymax for i,j in C.items()]
    H_min=[]
    H_max=[]
    sf = SIGMA[len_S-1]
    xH =[]
    for ci in x:
        ind_si = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==ci ]   
        H=[]
        for g_i in ind_si:
            G = istr2fsa(S[g_i],len_X,len_S)
            G1 = G.deletevent(sf)
            H.append(entropy(G1)/entropy(G))
        H_min.append(min(H))
        H_max.append(max(H))
    
    
    color2=(1,0, 0)    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)
    #l2=plt.step(xv,yv,".",linewidth=1.25,color=(0.5,0.5,0.5,0.5))
    plt.plot(x,y,"bo", color=color)    
    l1 = plt.vlines(x, [0], y, color=color,linewidth=1.5)
    l2 = plt.plot(x, H_min,"--", color=color2,linewidth=1.5,markersize=10)   
    #l3 = plt.plot(x, H_max,"--", color=(0,0,1),linewidth=1.5)
    # Format x-axis and y-axis.
    plt.axis([0, max_D, 0, max(y)*1.05])
    plt.xticks(range(0,max_D+1,5))
    plt.title(r"Entropy analysis of state complexity of $G_d$ for $||X||=%s$ and $||\Sigma||=%s$"%(len_X,len_S))
    plt.xlabel("Number of states of diagnoser\\\\ (worst theoretical case $2^{2\\times%s}=%s$)"%(len_X,2**(2*len_X)))
    plt.ylabel("Normalized frequency of instances")
    plt.legend([l2],["Minimum of entropy ratio"])
    f_path='../document/figures/Fig-entropy-%sn-%sk.pdf'%(len_X, len_S)
    plt.savefig(f_path, dpi=None, format='pdf')
    print "loaded Figure-entropy-%sn-%sk.pdf"%(len_X, len_S)
    if viz:
        command = 'evince'+' '+f_path+'&'
        retcode = check_call(command, shell=True, stderr=0)



def maximal_instances_D(n,k):        
    sf = SIGMA[k-1]    
    R = load('../instances/result-'+str(n)+'n-'+str(k)+'k.exp')
    print "result of experiment loaded"
    S = load('../instances/Gexhaustive-'+str(n)+'n-'+str(k)+'k.exp')
    #print len(S)
    print "instances loaded"
    max_D = R['max_D']
    ind_max_D = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==max_D]
    num_max_D = len(ind_max_D)
    Gmax = []
    for i in ind_max_D: 
        G = istr2fsa(S[i],n,k)
        G=G.setpar(Sigobs=SIGOBS[k])
        Gmax.append(G)
    save(Gmax,'../instances/GmaxDiag-'+str(n)+'n-'+str(k)+'k.exp')
    print "Gmax saved on hard disk"
    return

#instances=[(3,3)]#,(3,3),(3,4),(4,2),(4,3)]
#for i in instances:
#    len_X, len_S = i
#    analize_diagnosability(COLOR[len_X], len_X, len_S, viz =True)
    


def load_partial_results():
    len_V = []
    len_D = []
    ind_G = []
    test=[]
    res_db = shelve.open(RESULTFILE)
    for i in xrange(1,4):
        part = '%sn-%sk-%su-%sp'%(4,3,1,i)
        RP=res_db[part]          
        len_D += RP['len_D']
        len_V += RP['len_V']
        ind_G += RP['len_V']
        test += RP['test']
        instances = RP['instances']

    result_exp = {'n': 4,
                  'k': 3,
                  'instances': instances,
                  'valid': len(len_V),
                  'len_D': len_D,
                  'len_V': len_V,  
                  'ind_G': ind_G,  
                  'test' : test,
                  }
    fin = '%sn-%sk-%su-fp'%(4,3,1)
    res_db[fin]= result_exp
    res_db.close()
    print "saved file"
    
    


def generate_graphics_generic(color, instances, viz =False):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)
    for len_X, len_S in instances:    
         R = load('../instances/result-'+str(len_X)+'n-'+str(len_S)+'k.exp')
         print "result loaded for experiment -%sX-%sS"%(len_X,len_S)
         max_D = R['max_D']
         max_V = R['max_V']
         ind_max_D = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==max_D]
         ind_max_V = [R['ind_G'][j] for j, vm in enumerate(R['len_D']) if vm==max_V]
         num_max_D = len(ind_max_D)
         num_max_V = len(ind_max_V)
         num_diag = len([i for i in R['test'] if i==True])        
             
    
         H = Counter(R['len_D'])
         x = [i for i,j in H.items()]
         y = [j for i,j in H.items()]
    

         plt.plot(x,y,"bo", color=color)    
         l1=plt.vlines(x, [0], y, color=color,linewidth=1.5)
        # Format x-axis and y-axis.
         plt.axis([0, max(x)*1.05, 0, max(y)*1.05])
         plt.xticks(range(0,max_D+1,5))
         #plt.title(r"Frequency of state complexity of $G_d$ for $||X||=%s$ and $||\Sigma||=%s$"%(len_X,len_S))
         plt.xlabel("Number of states of diagnoser\\\\ (worst theoretical case $2^{2\\times%s}=%s$)"%(len_X,2**(2*len_X)))
         plt.ylabel("Valid instances")
        #plt.legend([l1],["$G_d$"])
         f_path='../document/figures/Fig-multiple-%sn-%sk.pdf'%(len_X, len_S)
         plt.savefig(f_path, dpi=None, format='pdf')
    print "loaded Figure-freq-%sn-%sk.pdf"%(len_X, len_S)
    if viz:
        command = 'evince'+' '+f_path+'&'
        retcode = check_call(command, shell=True, stderr=0)



    



#
#instances=[(4,3,2)]
#for i in instances:
#    print i
#    #calculate_diagver(*i)
#
#


#ins='%sn-%sk-%su-fp'%(4,3,1)
#    
#R=res_db[ins]
#max_D = np.max(R['len_D'])
#max_V = np.max(R['len_V'])
#ind_max_D = [R['ind_G'][i] for i, dm in enumerate(R['len_D']) if dm==max_D]
#ind_max_V = [R['ind_G'][j] for j, vm in enumerate(R['len_D']) if vm==max_V]
#num_max_D = len(ind_max_D)
#num_max_V = len(ind_max_V)
#num_diag = len([i for i in R['test'] if i==True])        
#Dm=np.mean(R['len_D'])
#Vm=np.mean(R['len_V'])
#Sm=np.std(R['len_D'])
#Sv=np.std(R['len_V'])
#print 'valid', len(R['ind_G'])
#print 'diagnosable', num_diag
#print 'max D', max_D
#print 'max V', max_V
#print 'mean D',  np.mean(R['len_D'])
#print 'mean V',  np.mean(R['len_V'])
#print 'maximum worst case D', len(ind_max_D)
#print 'maximum worst case V', len(ind_max_V)
#print "$(%s,%s,%s)$ & %s & %s &  %s & %.1f & %.1f & %s &   %s & %.1f & %.1f& %s \\\\ "%(len_X, len_S,num_unobs+1, R['instances'], len(R['ind_G']), max_D, Dm, Sm ,2**(2*len_X), max_V, Vm, Sv, 2*len_X**2)
#
#
#    




